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We investigate gas accretion flow onto a circumplanetary disk from a protoplanetary disk in detail 
, by using high-resolution three-dimensional nested-grid hydrodynamic simulations, in order to provide 

a basis of formation processes of satellites around giant planets. Based on detailed analyses of gas 
accretion flow, we find that most of gas accretion onto circumplanetary disks occurs nearly vertically 
toward the disk surface from high altitude, which generates a shock surface at several scale heights of 
the circumplanetary disk. The gas that has passed through the shock surface moves inward because 
its specific angular momentum is smaller than that of the local Keplerian rotation, while gas near 
\ the midplane in the protoplanetary disk cannot accrete to the circumplanetary disk. Gas near the 

O |. midplane within the planet's Hill sphere spirals outward and escapes from the Hill sphere through the 

two Lagrangian points Li and L2. We also analyze fluxes of accreting mass and angular momentum 
in detail and find that the distributions of the fluxes onto the disk surface are well described by 
power-law functions and that a large fraction of gas accretion occurs at the outer region of the disk, 
i.e., at about 0.1 times the Hill radius. The nature of power-law functions indicates that, other than 
the outer edge, there is no specific radius where gas accretion is concentrated. These source functions 
of mass and angular momentum in the circumplanetary disk would provide us with useful constraints 
on the structure and evolution of the circumplanetary disk, which is important for satellite formation. 
Subject headings: hydrodynamics — methods: numerical — planets and satellites: formation — pro- 
■ toplanetary disks — shock waves 
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1. INTRODUCTION 



Satellite systems around the giant planets in our solar system are commonly seen. They are thought to have formed 
in circumplanetary disks, which are believed to have existed around giant planets during their gas capturing growing 
stage. 

In earlier works, formation process of satellite systems have been considered based on a minimum mass subnebula 
(MMSN) model, in which satellites form from a disk that contains suffic ient solid mass with solar composition for 
reproducing t he current satel lite systems (e.g.. iLunine fe Stevenson! ()1982[ )). as an analog of the minimum mass solar 
nebula model (|Havashi 11198 If) . However, it was suggested tha t the MMSN model has difficulty in reproducing current 
satellite systems around Jupiter and Saturn ([Canup fc Wardll2002D . One of severe problems is that the model leads to 
much higher temperature than that of H2O ice sublimation at the current regular satellite region, which means that 
ice, which is the main component of the satellites, cannot be used as building material of the satellites. 

In order to overcome the dif ficulties of the MMSN- type models which assumes a closed and static disk, alternative 
models have been developed. iCanup fc Ward! ()2002j ) proposed a model in which an accretion disk with continuous 
supply of gas and solid is considered as a proto-satelli te disk. This model is based on results o f hydrodynamic simu- 
lations of gas capturing process of giant planets (e.g.. lLubow et al1ll999l ; iD'Angelo et al.ll200"2|) . which demonstrated 
that gas accretion from the protoplanetary disk toward the parent planets occurs through a circumplanetary disk. 
In this model, the surface density and temperature of the circumplanetary disk is kept lower than assumed in the 
MMSN model, and thus H2O ice can be used as a solid building material of satellites. Based on comparison among 
time scales of different processes, they concluded that formation of the Galilean satellites can be best explained if 
the circumplanetary disk had one order of magnitude lower gas surface density than the MMS N disk, with s low g as 
accretion rate corre sponding to Jupiter's growth time scale longer than a few xlO 6 yr (see also lSasaki et al.l (|2010ft ). 
On the other hand. iMosaueira fc Estradal (|2003T > proposed another disk model which consists of two components, i.e., 
inner MMSN-type massive disk and outer low-density extended disk. The model reproduced, for example, three inner 
Galilean satellites as well as only partially differentiated Callisto. Although these models seem to reproduce the current 
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satellite systems of the giant planets in our solar system, they needed to assume parameters for the disk structure such 
as surface density profile as a basis of physical processes of satellite formation. 

The structure of a circumplanctary disk is closely related to the gas accretion process of giant planets. Gas 
flow around proto-giant-plancts in protoplanctary disks has been studied using hydrodynamic simulatio ns. Earlier 
two-dimensional s i mulations showed that the gas i n the Hill sphere rotates in the prograde direction (|Mikil Il982l ; 
iSekiva et al.lll987t iKorvcansky fc Papaloizou Ill996f) and a pair of spiral shocks stands in the circumplanetary disks 
(e.g.. lKlevHl999HLubow et alJll999HTanigawa fc Watanabj|2002HD'Angelo et al.ll2002t) . However, the scale height of 
the proto-planetary disk for Jupiter-sized planets is comparable to the Hill radius, thus simulations with two-dimension 
approximation cannot capture the feature of the accretion flow. Recent three-dimensional hydrodynamic simulations 
reve aled that the two-dimensional picture for circumplanetary disks is not appropriate for the flow in the Hill sphere 
fe.g-.lD'Angelo et al.l2003t iBate et alJ2"00l iMachida et 1020081 iPaardekooper fc Mellernall2008t lAvliffe fc Batell2009bl : 
iCoradini et al.ll2010| ). One of the most important features newly found in three-dimensional calculations is that, in the 
region of circumplanctary disks, the gas accretes nearly vertically downward toward the midplanc from high altitude. 

Although some studies using hydrodynamic simulations were able to obtain the structures of circumplanctary disks, 
the direct use of results of hydrodynamic simulations for the structure, such as surface density, is problematic. One 
of the reasons is that disk-like shear-dominant flow is susceptible to numerical viscosity, which is intrinsic nature of 
numerical hydrodynamic calculation, no matter if SPH method or mesh based method is used, thus it is difficult 
to obtain reliable structure of the disk directly from such simulations. In particular, orbital radius of the rotating 
gas in disks basically changes only slightly with time, thus numerical error tends to accumulate easily in simulation 
of long-term evolution. Also timescale required for low viscous disks such as circumplanetary disks to reach steady 
state is much longer than the typical dynamical time of the fluid. While high-resolution hydrodynamic simulations 
is needed to resolve the structure of circumplanetary disks, their long-term evolution is difficult to follow with such 
time-consuming simulation. In addition, physical (non-numerical) viscosity in the disk is not well understood, and we 
need to assume specific viscosity models that are hard to justify. 

As an alternative approach, in the present work, we examine gas accretion flow onto circumplanctary disks from 
proto-planetary disks, in order to determine gas accretion rate as a function of distance from parent planets. Unlike 
the structure in a rotation disk, the accretion flow onto circumplanetary disks is not susceptible to numerical viscosity. 
Also, because the circumplanetary disks are located at the downstream of super-sonic accretion flow, the accretion 
flow itself is hardly affected by the circumplanetary disk structure, which depends on poorly-known effective viscosity. 
In Sj2j we describe our settings of hydrodynamic simulations. In Sj3j results and analyses of the simulations are shown. 
We discuss implication of our analyses of the accretion flow in 2) We summarize our results in §5\ 



We consider a situation in which a growing giant planet embedded in a protoplanctary disk has induced the nucleated 
instability and the gas of the disk accretes dynamically onto the planet. We take local Cartesian coordinates rotating 
with the planet, a;- axis corresponding to the radial direction, y-axis the revolving direction of the planet, and z-axis 
normal to the disk midplane. The planet is located at the origin. 

We adopt the local approximation, in which tidal potential is linearized and curvature of the protoplanctary disk is 
neglected. This approximation is valid as long as the Hill radius is much smaller than the orbital radius of the planet. 
The orbit of the planet is assumed to be fixed circular and co-planar with the disk midplane. The gas is assumed to 
be inviscid and isothermal. Magnetic field and self-gravity of the gas are neglected. 

We employ equation of continuity, equation of motion for compressive inviscid gas, and equation of state for isother- 
mal gas to simulate the gas motion. In order to understand physics clearly, we use these equations in a non-dimensional 
form. We normalize time by the inverse of the Keplerian angular velocity ilj^ 1 = (GM*/a 3 ) -1 / 2 , length by scale height 

h = cSl^ 1 , and mass by unperturbed surface density of the protoplanetary disk gas So times ft 2 , where G is the 
gravitational constant, M* is the solar mass, a is the semi-major axis of the planet, and c is the sound speed of the 
disk gas on the planet orbit. The sound speed c is unity in our normalization. Normalized quantities are denoted with 
tilde on the variables in this paper. The normalized equations can be written as 



2. METHODS 



2.1. Settings and Basic Equations 
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*p = (6) 

In the above f = \ r\ is distance from the origin, fn is normalized Hill radius, fjj = r^/h = (M p /3M ii ) 1 ' 3 (a/h), M p 
is the planet mass, and e~ = (0, 0, 1) is a unit vector in the z-direction. The third term of the right hand side in 
Eq. ([!]) is added so that <f> = at the two Lagrangian points Li and L2 (i.e., (x,y,z) = (±fH,0,0)). The normalized 
quantities can be written, for example, as v = v/(hflK) = v/c, p = p / (Y<q / {\f2~Kh)) , P = -P/(c 2 (So/(v / 27r/i))), and 
$ = $/c 2 f or velocity, den s ity, p ressure, and potential, respectively. Full equation set before the normalization can 
be found in iMachida et al.l ()2008| ). In the normalized system, rn is the only characteristic physical parameter of the 
system, if we do not consider physical radius of the planet. In this paper, we focus on the case with fn = 1, which 
roughly corresponds to Jupiter mass at 5AU. 



2.2. Numerical Modeling 

We employ a three-dimensional nested-grid hydrodynamic simulation code (e.g.. IMachida et aLll2 005, 200 6)), which 



was originally developed to explore star formation process by collapse of molecular cloud core (jMatsumoto k, Hanawal 
12001 . Size of the whole computational domain (L x , L y , L z ) is (24,24,6). The domain has a symmetry about 2 = 
plane, and the covered region in simulation is x = [—12, 12], y = [—12, 12], and z = [0,6]. We set 11 levels for the 
nested-grid system, and the numbers of grids in each level are (n x ,n y ,n z ) = (64,64, 16). The level of the nested grid 
is denoted by I and / = 1 is the largest grid level. Increment of I by one reduces the size of the computational domain 
to half in all directions, keeping its center at the planet. The finest grid size is thus L x /64/2 11 ~ 1 = 0.000366, which 
corresponds to about 1/4 of the present Jovian radius at 5.2 AU. 

We adopt two types of artificial gravitational weakening, which allows us to avoid non-physical crash of the com- 
putation without essentially changing results. One is the widely-used weakening for region adjacent to the planet to 
avoid singularity of the planet gravity. The modified potential is given by 



3f| 



)l/2 

1/ 



(7) 



where r sm is so-called smoothing length and we set f sm = 0.00073, which corresponds to two grid size of the finest 
grid level (I = 11). The other gravitational weakening is also applied for the z-component of the tidal force in the 
high-z region. In this region, tidal force in the z-direction (toward the mid-plane) is very strong, and causes too large 
density gradient to be described by the grid size we set in the code. Thus we connect two parabolic curves smoothly 
at z = L* 



as 
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where we set L* = 4. Initial density profile in the z-direction is modified accordingly as 
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(9) 



for hydrostatic equilibrium under the potential given by Eq.®. This weakening does not cause any significant effect 
for the region where we need to observe because the mass above L* is negligible (0.0063 % of the total mass for the 
hydrostatic structure). 

As boundary conditions, we set the flow at i = ±L x /2 to be the unperturbed one, and the mirror condition is applied 
at z = and L z . As for the boundaries at y = ±L y /2 and i ^ (inflow region), we use mixed boundary conditions as 
follows. When t = 0, we set unperturbed condition, which is the same as the initial condition, and change it gradually 
to periodic boundary condition until t = 100 based on linear interpolation with respect to time. When t > 100, 
periodic boundary condition is applied to all region at y = ±L y /2. For the initial condition, we adopt unperturbed 
density and velocity; v x = v z = 0, v v = — (3/2)x for veloci ty and dp/dx = dp/dy = 0, p(z) = {2-k)~ 1 / 2 cxp(— z 2 /2) for 
density, which is identical to that of IMachida et al.l ()2008l) . 

In order to obtain high-resolution results efficiently, calculations are started with only the largest grid level and the 
number of levels is increased with time. The second largest level (/ = 2) starts at t = 20, I = 3 at t = 50, and the 
higher levels (/ > 4) are started in order from t — 150 by considering relaxation degree in each level. In this paper, we 
will show a snapshot of the flow of t = 160.7, when the accretion flow is nearly in the steady state. 
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We set sink cells around the origin in order to see the pure effect of gas accretion flow without the effec t of the 
planet body, as well as to mimic gas accretion phase ont o the planet that follows the nucleated instability (e.g.. lMizunol 
119801: iBodenheimer fc Pollaclj|1986t llkoma et al~ll2000D . In the sink cells, gas is removed at a rate that corresponds 
to p/p = 10~ 4 . We set f sink = f sm . 

3. RESULTS 

3.1. Path to Circumplanetary Disks 

First, we examine the overall structure of accreting gas flow onto a circumplanetary disk. Figure [T] shows streamlines 
from outside of the Hill sphere (i.e., from a protoplanetary disk) toward the planet. The starting point of a streamline 
is defined by (xo,yo,zo)- Four panels show streamlines starting from four different height; So = 0.0, 0.5, 1.0, 1.5. 
Based on the destination of streamlines, we divide the flow into three regions: passing region where gas approaches 
the planet and passes by it without making U-turn, horseshoe region where gas crosses the planet orbit with U-turn 
during the encounter, and accreting region where gas is accreted onto the circumplanetary disk. We find that there are 
no accreting streamlines on the midplanc (z = 0), while some streamlines starting from off-midplane region accretes 
to the circumplanetary disk. This implies that the gas near the midplanc in proto-planctary disks is harder to accrete 
to circumplanetary disks and planets. This is confirmed by Figure [2j which shows streamlines starting from three 
different radii R = ^/x 2 + y 2 = 0.05, 0.1, 0.2, that is, (x ,yo,z ) = (0.05, 0, 0), (0.1, 0, 0), (0.2, 0, 0). We can see that 
gas on the midplanc in 0.2 < R < 1 spirals outward and escapes from the Hill sphere within a short timescale through 
one of the two Lagrangian points (Li or L2), and that the outward radial velocity decreases with decreasing distance 
from the planet (see also Fig. [7]). 

The vertical structure of the flow is more clearly demonstrated in Figure [3] In this figure, initial positions of 
streamlines on the x-z plane at y = L y /2 are classified into three regions by their destinations. We define the passing 
region where streamlines reach the boundary at y = —L y /2 with x > 0, horseshoe region where streamlines reach the 

boundary at y = L y /2 with i < 0, and accreting region where streamlines end up within a sphere with radius ft,; 
we set fb — 0.2. We can clearly see that there is no accretion band in the midplanc, while off-midplane site has an 
accretion band with a significant width. Note that the classification of the three regions does not depend on fb as long 
as 0.1 < f b < Ly/2. This implies that, even in the Hill sphere, the gas in the region f > 0.1 stays there temporarily 
and is not necessarily captured by the planet. We will discuss the flow in this region in £)4.1I 

3.2. Disk Structure and Gas Motion 

Figure [4] shows azimuthally averaged density, velocity, and specific angular momentum of the gas with three different 
scales. We first recognize that the gas clearly forms a disk-like structure; the density distribution is concentrated near 
the midplane, almost Kcplcrian rotation is realized in the high density region (judging from the contour lines of specific 
angular momentum), and the flow velocity in the R-z direction of the region is very low. The density profile in the 

z-direction can be roughly described by hydrostatic equilibrium in the region where z < bh v and R < 0.2, where h p is 
scale height of the circumplanetary disk defined by 

h p = ^f¥J(3^)o,R^ 2 . (10) 

From this dependence a flare-up disk is expected, which is seen indeed in Fig.|4j Above the disk surface, large downward 
velocity, which is almost free fall velocity, is observed; this shows that gas is indeed accreted directly onto the disk 
surface. The accreting gas forms a shock surface, which stands at z ~ 5h p . By analyzing streamlines, we confirmed 
that these gas elements are actually originated from off-midplane gas (mostly from z > 0.5) in the protoplanetary 
disk. Also, the contour lines above the shock surface are aligned with the velocity vectors, which means that specific 
angular momentum does not change very much before the gas elements hit the disk surface in this inner region, where 
three-body effect is weak enough to be neglected. 

Figure [5] shows radial mass flux pv r (where v r = v ■ r/f), as a function of the azimuth and elevation angles at 
spheres with four different radii f = 1.0, 0.3, 0.1, 0.03. This figure clearly shows that the accretion manner strongly 
depends on both these angles and is not spherically symmetric at all. For example, the flux at the f = 1.0 sphere 
shows that the mass flux can be both inward and outward directions near the midplanc (where \9\ < 40°, where 9 
is elevation angle) depending on azimuth angle while it is always inward (i.e., negative flux) at high \9\. The two 
maxima near </> = and 180 on the midplane {9 = 0) in the case of f = 1.0 corresponds to the outflow from the 
Hill sphere through the two Lagrangian points Li and L2 shown in Fig. [2] The flux on the other three spheres also 
shows that there are two positive maxima and negative minima near the midplane, which implies the formation of a 
two-arm spiral structure in the circumplanetary disk. The range of the elevation angle 9 that corresponds to outward 
flux shrinks with decreasing f, which indicates the change of disk aspect ratio with f . 

Although Figure [5] shows that the gas at high elevation angle is always falling radially inward, it is difficult to judge 
the direction of net mass flux in the disk (corresponding to low \9\) from these plots. In order to examine the net mass 
flux, next we calculate azimuthally integrated mass flux at spheres with several radii as a function of 9 (Figure EJa)) : 



F r (r,6) = / {pv r y 2 cos 8dcj). 
Jo 



(11) 
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We find that the direction of the net radial flux near the midplanc (\6\ < 20° in the case of f = 1.0), where the gas 
moves both inward and outward through the sphere in the region depending on the longitude, is actually outward, 
while it moves inward for higher elevation angle (16*1 > 20° in the case of f = 1.0), which is obvious in Fig. [5] This 
basic profile does not change with f, but where f < 0.1, the profile becomes sharper, that is, the region where gas 
moves outward is limited in narrower elevation angle, which corresponds to the dependence of disk aspect ratio on 
radius. This clear outward motion around the midplane for wide range of f shows another evidence of the unbound 
state of the gas where f > 0.1, which is already shown in the above (Fig. In the cases with f = 0.1 and 0.03, there 
is a narrow band where the mass flux is inward (negative) with peak at 6 ~ ±25° and ±15°, respectively. Figure^b) 
shows azimuthally averaged radial velocity: 

s /- m_ J 2 "p(r,MK(r,fl,</>)# . . 

Jo p{r,e,4>)d4> 

This quantitatively shows that the velocity of the gas elements accreting toward the disk nearly vertically at high 
elevation angle (\0\ > 40° in the case with f = 0.1) is nearly free-fall velocity (~ \/6f : j [ /f), which is faster than the 
sound speed (= 1). Near the midplanc, radial velocity is very small, which corresponds to the hydrostatic disk region 
described above (Fig. [4}. Between the two regions, there is a transition region, which is under the shock standing at 
the disk surface, and the gas in the region is not in dynamically equilibrium state as a part of the disk. The range of 
9 of the transition region corresponds to the narrow band with inward mass flux described above. Thus, this is a kind 
of layered accretion and its mechanism will be discussed in § 13.31 
Figure [7] shows radial velocity in the midplane normalized by local Keplerian velocity defined by 

v.(fi)-M^M=a. (is) 

We find that radial velocity is very small but takes on positive values for a wide range of radii (0.005 < R < 0.2). 
Significant positive values at R > 0.2 correspond to the outward motion which is already shown in Fig. [2j Non- 
negligible positive values of Vr at R < 0.005 may reflect non-steady state small-scale structure in the region, or 
artifact of averaging operation over discrete grids whose size may not be sufficiently small as compared to the distance 
from the origin. Formation of the circumplanetary disk with nearly Keplerian motion is also demonstrated in Fig. [8j 
which shows azimuthally averaged specific angular momentum at the midplane. We clearly see that the rotation 
velocity is nearly Keplerian for a wide range of radii (R < 0.1). Keplerian rotation is realized when the ratio of 
thermal energy of the gas to potential energy of the planet is much smaller than unity (or pressure force is much 
weaker than gravity force). Eq. (|10j) indicates that h v /R, which is square root of the ratio, becomes smaller with 
decreasing i?, which is consistent with our result. 
Figure [5] shows the plots of azimuthally averaged density at the midplane and surface density. We can see that the 

profiles at R > 0.02 can be fitted roughly by power-law functions as p oc R~ 3 and S oc i?~ 3 / 2 , respectively. In an 
equilibrium state, we have p = E/(-\/27r/Zp). The above two functions imply that this relationship is approximately 
satisfied and the hydrostatic equilibrium is realized. Deviation from the power-law functions for the inner region may 
be due to insufficient computation time as compared to the time required to reach a steady state 1 , or due to the effect 
of sink cells around the center. Here we emphasize again that the distributions of mass and angular momentum of 
accreting gas that we present in § I3.3 | are more impo rtant and useful quantities less affected by numerical procedures, 
as we mentioned in § [TJ Note that IMachidal ([20091 ) introduced an idea of centrifugal barrier to explain the peak of 
surface density. However, we do not find any physical reason for the existence of such a barrier because, as we will 
show in £13-31 distributions of accreting mass and angular momentum onto the circumplanctary disk are well described 
by power-law functions, which do not have typical lengths, such as centrifugal radius. 

Figures [TU] and [11] show variation of physical quantities along two streamlines starting from different heights that 
correspond to the passing and accreting regions. Overall variation is shown in Fig. I10[ while close-up view near the 
shock front is shown in Fig. 1111 The key quantity in this figure is the Bernoulli integral B given by 

B = i| v\ 2 +\ogp + $, (14) 

where we included potential energy. This quantity is conserved along each streamline except shock surfaces. 
First we examine a streamline in the midplane plane (Figs. [TUJa) andQTJa)). The starting point of the streamline 

is (xo,yo,Zo) = (2.37,2^/2,0), which is in the passing region (red region in Figs. [1] and [3]); x keeps positive values 
not far from 2, and y monotonically decreases with increasing d, where d is the distance from the starting points of 
streamlines. With approaching to the planet (corresponding to d < 11.6), the gas element climbs the tidal potential 

1 Time required to reach steady state can be estimated as (7rii 2 S/Af) r= o.o2 ~ 30, which is longer than the duration with full level 
calculation. Here we used M = M B (R) (see Fig. 1 1411 . 
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slope decreasing its kinetic energy, and density slightly decreases. Although these quantities change with d, Bernoulli 
integral B, which is the sum of the three quantities, is constant until d ~ 11.6, which shows there are no shocks along 
the path until this point. However, there is a shock surface at d ~ 11.6 where kinetic energy decreases and density 
increases discontinuously, and B decreases as a result of shock dissipation. This shock surface stands from the Hill 
sphere toward both inside and outside of the planet orbit, forming spiral density waves around the central star. After 
passage of the shock, B remains nearly constant, although other quantities change with increasing d. This means that 
the gas element on the streamline underwent only one shock with modest strength. 

Next we examine an off-midplane streamline (Figs. HOf b) and lllf b)). The starting point is (io,yoj^o) = 
(2.37,1/^/2,1), which is just above the previous starting point with the same horizontal position, and this corre- 
sponds to the accreting region (blue region in Figs. [T]and|3]). The quantities change in a similar way to the case in 
the midplane before d ~ 12.0, where the gas element hits the shock surface forming around the Hill sphere. We have 
non-zero z values in this case but it is kept almost constant until it encounters the shock, which shows that the flow 
before the shock is laminar. When the gas element passes through the shock surface at d ~ 12.0, basic behaviour is 
the same as the zq = case, but the gas element changes the direction slightly upward (positive z direction) at the 
shock (see the line showing the change in z in Fig. [TTT b)). This is because the shock surface forms a bow shock and is 
curved upward, thus the gas element of off-midplane, which passed an oblique shock, gains upward momentum. After 
that, the gas element falls steeply toward the mid-plane and hits the surface of the circumplanetary disk, forming 
another shock at d ~ 13.9. The point of fall is at f ~ 0.1, which is sufficiently close to the planet to make use of its 
gravity for acceleration. Thus Mach number of the gas when it hits the disk surface is much larger than unity, and 
strong shock dissipation is caused. Owing to this strong dissipation, the gas becomes captured within the Hill sphere 
and merged as a part of disk. Note that the gradual decrease of B before d ~ 13.9 arises presumably from collisions 
between falling gas elements on the way to the surface, especially between two gas elements originated from x > and 
x < in the protoplanctary disk (see also Fig. [T2l and the corresponding description in the text). After the passage 
of the shock forming above the circumplanetary disk (d > 14) , f becomes almost constant and azimuth angle changes 
with a constant rate, which means that the gas is now rotating around the planet, i.e., the gas element has become a 
part of the circumplanetary disk. Note that specific angular momentum is not constant along the streamlines at all, 
although it is sometimes assumed to be a conserved quantity. 

Now we consider why the gas in the midplane cannot accrete to the circumplanetary disk. As described above, 
effective shock dissipation is required for gas elements to become accreted onto the circumplanetary disk, and the 
planet gravity is essential for the effective acceleration. When gas element consumes potential energy without significant 
increase of gas density, kinetic energy of the gas increases effectively (see Eq. ([14)) ) and the gas element can have a 
high Mach number, which leads to strong shock dissipation. However, in order for the gas in the midplane to consume 
planet gravitational energy, the gas element has to pass through the high density region of the circumplanetary disk. 
The increase of density absorbs potential energy and prevents gas from effective acceleration. Gas near the midplane is 
thus difficult to form strong shocks and difficult to dissipate energy effectively, which prevents it from being captured 
by the planet gravity. Therefore, the gas in the midplane is difficult to be accreted into the circumplanetary disk. 

One may think that passing many weak shocks during rotating around the planet in the circumplanetary disk would 
contribute to inward migration of the rotating gas. However, the decrease in B is a third-order small quantity with 
respect to (Ai — 1) (where M. is Mach number), so such weak shocks are not effective for energy dissipation and 
inward migration. Note also that acceleration by the tidal force in the z-direction alone is not sufficiently strong 
to form strong shocks. This is because the gas is in hydrostatic equilibrium supported by thermal pressure, so that 
the resultant accelerated velocity should not be much larger than the sound speed, being able to produce only weak 
shocks. We therefore conclude that the strong shock produced by direct infall onto the surface of the inner part of the 
circumplanetary disk and associated energy dissipation is essential for the accretion of the gas from protoplanetary 
disk onto the circumplanetary disk. 

3.3. Distribution of Accreting Mass and Angular Momentum onto the Disk 

Next, we examine accretion rate of mass and angular momentum as a function of distance from the planet. First we 
define the surface at which physical quantities immediately prior to the accretion onto the surface of a circumplanetary 
disk is measured. The surface is axisymmetric about the z-axis and the height of the surface from the midplane is 
denoted by z s (R). We define the height of the surface by a linear function of the local scale height of the circumplanetary 
disk as 

z s (R) = f h h p + 2 s0 , (15) 

where //, and z s o are constants. We set fh = 6 and z s q = 0.004 so that the surface z s is above but not far from the 
shock surface of the circumplanetary disk in the region with R < 0.2. 
Next we define mass flux onto the circumplanetary disk by 

/,(^^)= >' t, '. n> , (16) 

COS U s 

where p a and v s are the density and velocity vector at the surface defined by Eq. (|15p . n s is an upward unit vector 
normal to the surface, 8 S = dz s /dR, and the negative sign is added so that f s becomes positive when the gas passes 
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the surface downward. Note that cos 9 S appears in the denominator in Eq. (|16|) because / s is defined as a flux at the 
surface per unit area of the x—y plane, not a flux per unit area of the surface. 

Figure [12] shows several quantities on the surface z = z s projected on the x-y plane. We find that the density 
of the accreting gas at the surface (p s ) has peaks at (x,y) ~ (±0.1, =p0.08), while the mass flux / s has peaks at 
(x,y) ~ (±0.04, =p0.04), which are located closer to the planet than the density maxima. This is because accreting 
velocity, which is almost free-fall velocity, is faster at the inner region. These two peaks of the two-arm structure 
correspond to the places where gas elements coming from interior (x < 0) and exterior (x > 0) to the planet orbit 
collide with each other. As for the z-component of specific angular momentum around the planet measured in the 
inertial frame j z>s = ( r X t>) a tz=z 3 • e 2 + R 2 (where the R 2 term arises from the rotation of the coordinate system), its 
value is positive (i.e., progradc) in the whole region. The distribution is roughly axisymmetric with some elongation 
in the x-direction within R < 0.1, and the value increases with increasing radius. As for the angular momentum flux, 
the above axisymmetric property breaks and has peaks at (x,y) ~ (±0.12, =f0.08), which simply corresponds to the 
two-arm structure observed in density or mass flux. Note that the negative values observed in / s in the outer region 
(R > 0.1) are due to the fact that the gas passes through the z = z s surface from beneath. This is because, in this 
region (R > 0.1), the ratio of the scale height to the radial distance of the circum-planetary disk is not much smaller 

than unity anymore and z s is comparable to R accordingly. This leads to negative / s for a gas element whose horizontal 
velocity is comparable or larger than vertical velocity even when the vertical velocity is negative (downward). The 
surface of the disk in this region is hard to define, thus the way to measure mass flux described above may not be 
appropriate in this region. 

For a better understanding of the distribution of accreting gas, we next show azimuthally-avcragcd quantities at the 
surface defined by Eq. (fl~5j) . First we examine mass accretion rate. Figure [T3] shows the plots of azimuthally-averaged 
mass flux onto the disk, defined by 

1 f*« - _ 

f B (R) = — f a (R,<f>)d4>. (17) 

Figure [T4l shows the plots of the cumulative mass accretion rate, i.e., the mass accreted per unit time through the part 
of the surface z s where horizontal distance from the planet is smaller than R, i.e., 

rR p2lT 

M S (R)= / f s {R,4>)R'd(pdR'. (18) 
Jo Jo 

We can see from Fig. [T3] that f s (R) is almost constant from 7? ~ 0.1 all the way to the very center (R ~ 0.001). 

Cumulative mass accretion rate M s is thus roughly proportional to R 2 , which is a power- law function and thus does 
not have any typical length except for the outer end of the downward accretion on the disk surface (at R ~ 0.1 in this 
case). This distribution shows that the main contribution of mass accretion comes from the outer region (R ~ 0.1) 
and there is no specific radius where gas accretion is concentrated somewhere in between. Note that dotted lines in 
Figures H31 and UM show quantities calculated by Eqs. ([XT]) and (fTgj) , respectively, integration being performed over 
only regions with / s > 0. Deviation from each solid line can be only seen at R > 0.15, where gas can pass through the 
z — z s surface even upward depending on azimuth angle (see Fig. [T2|) . This is the reason why M s (solid line) decreases 
with increasing radius at R > 0.2 although it is a cumulative quantity. The treatment of our analysis in this region is 
a difficult problem and will be discussed later f H4.ll) . 

The distribution of / s (or M s ) provides us with useful information but it is not sufficient for an understanding of 
the mass distribution in a circumplanetary disk. In general, the radial distance from the planet of the location where 
gas is first accreted on the disk is different from its final orbit after it settles in as a part of the Keplerian rotating 
disk. What connects the two is angular momentum of the gas. Figure [15] shows azimuthally-averaged specific angular 
momentum j z _ s at the z = z s surface given by 

/ 3z,sP s d(f> 

l-AR) = jjL jr« ' (19) 

/ p s d(j) 
Jo 

We can see that j z . s is always smaller than the specific angular momentum for Keplerian rotation jKcp = \J 3r^R. 
This means that the gas does not have enough angular momentum for Keplerian rotation at the radius, which would 
lead inward migration of the gas after the fall on the disk surface, while the gas near the midplane moves outward, 
as mentioned above (Fig. [7]). This explains the inward stream of the layer just under the shock of the disk surface 

(Fig. [S]). In addition, j zs is nearly proportional to R, which is steeper than jKep cx R 1 ' 2 , and thus the ratio of j ZjS to 
jKep decreases with decreasing R. This suggests that the gas accreted in the outer region of the disk (R ~ 0.1) tends 
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to keep the position when it rotates as a part of the circumplanctary disk, whereas the gas accreted in the inner region 
tends to move further inward in order to achieve balance between centrifugal and gravitational forces. 

We now have the mass flux and angular momentum in the accretion flow, which allows us to estimate the effective 
distribution of accreting gas elements, assuming their redistribution to radial distances where their specific angular 
momentum matches that of the local Kcplcrian rotation. Let i?Kep denote the radius where the gas with specific 
angular momentum j z s is rotating with the Keplerian velocity, i.e., 

RKe P (R,<l>) = £t ■ ( 2 °) 

Then the effective distribution of azimuthally-averaged mass flux after redistribution of the gas with angular momentum 
conservation can be written as 

f Kep (R)dR=—=r dR' R'd<f>f s (R', 4>)&{R' , 4>), (21) 
2ttR Jo Jo 



6(R' <j>) = { 1 iiR - dR / 2<R ^ P {R'^)<R + dR/2andf s >Q, ^ 



where 

A(R' rh) = J 

otherwise 

The corresponding cumulative mass accretion rate is 

M Kcp (R) = / 2irR'f Kep (R')dR'. (23) 
Jo 

The plots of /Kep and MKep are shown in Figures IT3l and fT4l We can see that /Kep is roughly proportional to R" 1 
and Mkcp is to R. By comparison with / s , we find that fuep is more center-concentrated distribution, but still outer 

region is a dominant source of gas accretion in the sense that d In fKep/d In R > —2. Since angular momentum of 
a gas element after the shock of the disk surface is not necessarily conserved until the gas reaches the radius where 

it rotates in the Keplerian velocity, the distributions of /kcp and Afxcp should be regarded as the case of the most 
center-concentrated limit. 

4. DISCUSSION 

4.1. Size of Circumplanetary Disks 

The size of a circumplanetary disk is closely related to the location of satellite formation. However, it is not easy 
to define the outer edge of the disk, and there have been several attempts. One natural way would be to define the 
disk edge based on its density distribution. But the density is monotonically and smoothly d ecreasing with increasin g 
radius toward the Hill sphere, so it is difficult to define the edge from density distribution. lAvliffe fc Bate! (|2009b[) . 
who performed three-dimensional hydrodynamic simulations, suggested a criterion for the disk edge based on angular 
momentum distribution. Their simulations show that specific angular momentum of the disk gas has a peak at a 
certain radial location, whereas specific angular momentum of a Keplerian disk increases monotonically. Their Figure 
2 suggests th at the turnover point is about 1/3 of the Hill radius, and they defined the disk edge by the radial location 
of this peak. iMartin fe Lubowl (|2011l ) examined periodic orbits of a particle around a planet under the influence of the 
gravitational force from a central star and the planet, and found that, as the size of the orbit is increased, the orbits 
start crossing with each other at r ~ 0.4rn. They inferred that this corresponds to the location of the disk's outer edge 
where tidal torque of central star's gravity becomes strong, and called it tidal truncation radius (r trunc ). The above 
value of r trU nc is in agreement with the point of turnover of specific a ngular momentum found in their two-dimensional 
SPH simulation, as well as in the result of lAvliffe fe Bate! (|2009bf ) . Our Figure [8] shows that the turnover point is 

R ~ 0.3rn, which is also roughly in agreement with rt rU nc-* 2 

We here consider an alternative criterion by examining Vr defined by Eq. (|13l) to define the position of the outer 
edge. We observed that azimuthally-averaged radial velocity in the midplane is positive (outward) in almost all the 
region, although the value is very small for R < 0.2 (Fig. [7]), and the outward velocity significantly increases at R > 0.2 
(see also Fig. Since gas at R > 0.2 moves outward and escapes from the Hill sphere quickly, the region R > 0.2 
can be regarded as outside of the circumplanetary disk. We thus define the disk edge as the radial location where Vr 
starts increasing significantly and has non-negligible positive value; R ~ 0.2 in our case sho wn in Fig. [7J Note tha t 
a two-dimensional SPH simulation shows that negative torque is exerted on the outer disk ([Martin fc Lubowl 120 111 ), 
which gives rise to inward flow and is contrary to our results. This might be due to the fact that two-dimensional 
simulations, which create clearer spiral structure than three-dimensional ones, tend to enhance torque density. In 

In this comparison, wc assume that specific angular momentum in these works were not measured in incrtial frame, but in the rotating 
frame. 
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particular, the disk thickness in the outer region (R > 0.2). at which outflow is observed, is very thick (thickness 
is comparable to radius). Thus, it would be unlikely t hat o ur three-dimensional calculation is significantly affected 
by the negative torque suggested bv lMartin fe Lubowl (|201 lh . The distribution of accreting angular momentum also 

seems to indicate a similar radius for the edge. As we see in Fig. 1151 j z , s and jKep have different dependence on R. If 

we fit each profile as a single power-law function, they cross each other at R ~ 0.3 in the present case. This roughly 
agrees with the location of the disk edge defined above based on the significant increase of outward velocity. The gas 
exterior to this radius has too large angular momentum to achieve Keplerian rotation, thus moves radially outward. 

Once angular momentum of accreting gas is obtained, one can calculate mean specific angular momentum of the 
accretion flow, which is sometimes used in calculating the the s o-called centrifugal radius f c = £ 2 /(3f^) to infer the 
disk size (e.g., Mosauei ra fc Estradal 120031 iWard fc Canupll2010( ). Mean specific angular momentum of the accretion 
flow within a radius R is given by 

i(R) = J^i, (24) 
M S (R) 

where 

pR p-2k: 

J S (R)= dR! R'd<t>f s {R'A)UR'A)- (25) 
Jo Jo 

We can see from Figure [T51 that t{R) increases nearly linearly with radius and levels around R ~ 0.2, where I ~ 0.2. 
This corresponds to f c ~ 0.013, which might give us a rough estimation of the disk size. However, it is practically 
difficult to determine the upper bound of the integral range with respect to R. which affects the value of £, because the 
circumplanctary disks is smoothly connected to the proto-planetary disks. Also, unlike the case of a particle, angular 
momentum of a gas element is not a conserved quantity along a flow in principle, because gas is continuum medium 
which can transfer angular momentum through waves. This also makes it difficult to determine i precisely, and thus 
f c as well. 

4.2. Picture of Gas Accretion Flow onto Circumplanctary Disks 

Figure [T71 shows a schematic picture of gas accretion flow onto a circumplanctary disk based on the results obtained 
in the present work. Gas accretion occurs mostly downward from high altitude with high incident angle (see Fig. 2]). 
The accreting gas is accelerated by the planet gravity to have almost free-fall velocity of the planet (Fig. [Bb) . The 
value of angular momentum of the accreting gas normalized by the local Keplerian angular momentum is lower in the 
inner region of the disk and higher in the outer region (Fig. I15|). The falling gas reaches the shock surface formed on 
the top of the circumplanetary disk. 

Gas near the midplane (especially where R < 0.1) is almost in Keplerian rotation and hydrostatic equilibrium in 
the z-dircction. In this region, radial velocity is very small and it is difficult for the gas to accrete inward through 
the disk midplane, as mentioned above. On the other hand, the gas at R > 0.2 in the midplane shows significant 
outflow, which eventually escapes from the Hill sphere (see Fig. [2]). Thus inflow from high altitude and outflow near 
the midplane co-exist, and they do not interfere with each other. This means that there is a circulation across the 
Hill sphere and fresh (protoplanctary) gas is always supplied in the R > 0.2 region. In addition, the two Lagrangian 
points Li and L2, which are often thought to be the most likely po ints of inflow, are ac tually the points of outflow, 
even in the gas accretion stage (Fig. [5]). This seems consistent with lKlahr fc Klevl (|2006l ). who also performed three- 
dimensional hydrodynamic simulation with no explicit physical viscosity, but with lower resolution and with radiation. 
They showed that accretion mainly occurs via the poles of the planet and no inflow along the equatorial plane, which 
is quite similar to ours. But they explained that the circulation in their simulation is driven by accretion heating, 
which we do n o t consi der. Thus it is not clear if the driving mechanisms of the circulation are the same. Note that 
lAvliffe fc~B ate (2009b) shows that the vast majority of the mass flows into the Hill sphere near the equator, which 
seems to be inconsistent with our result. There are several differences in setting between theirs and ours, so it is not 
easy to judge which factor causes the difference. Disk thickness would be one of the reasons for that, even though our 
disk thickness is not very different from theirs. Another possibility is gap formation. When a deep gap is formed, the 
contribution of the vertically accreting gas would become less significant, which reduces the difference. However, here 
we think that viscosity is the main reason for the difference. When there is explicit viscosity, the circumplanetary 
disk gas should transfer its angular momentum outward and most of gas would move inward accordingly, but this is 
not necessarily true in the inviscid limit. Actually, inviscid simulations bv lKlahr fc Kiev! (|2006| ) showed similar results 
with ours. Also, invisc id limit might not be bad because circumplanetary disks are likely to be MRI inactive in most 
cases (jFuiii et al.ll2011[ ). 

Since dust particles tend to settle down toward the midplane, gas accretion flow from high altitude is likely dust-poor 
gas, which diminishes dust-to-gas ratio in the circumplanctary disk. On the other hand, since gas in the outflow region 
(R > 0.2) rotates significantly slower than Keplerian velocity around the planet, dust particles would migrate inward 
quickly. Thus the outflow in the midplane is also likely dust-poor gas, which would be the source of solid material in the 
circumplanctary disk and would enhance dust-to-gas ratio in the disk. Dust-to-gas ratio in the circumplanetary disk 
is one of the most important factors for satellite formation processes (jCanup fc Wardll2006| ). thus these two filtering 
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effects would become import ant for satellite forma t ion. In addition, studies on size evolution of solid material in 
protoplanctary disks, such as iKobavashi fc Tanaka I (|2010l) . would also be important for the filtering effects. Further 
studies are needed on this issue. 

As we see in Figure |6l there is the layer just under the shock at the disk surface where gas moves inward. This 
layered accretion is explained by the fact that angular momentum of the accreting gas onto the disk surface is lower 
than that of the Keplerian rotation at the radius (Fig. I15j) . The inward stream in the layer might play an important 
role for the net mass accretion toward the planet in circumplanetary disks. This accretion picture should be explored 
more quantitatively by further high-resolution simulations, using a code with lower artificial viscosity such as the one 
with polar coordinates for numerical grids. 

4.3. Effects of Gap Formation 

One important issue to be addressed is the effect of gap, which is a lower density annulus region near the planet 
orbit in the protoplanetary disk. We observed the gap as a slightly lower density band formed around x = in our 
simulations. However, in the last stage of giant planet formation, a giant planet would become massive enough to 
create a deep gap, which would be able to truncate its growth. Since the deep gap is associated with steep density 
gradient at the edge, gap formation may affect accretion flow and circumplanetary disk formation. However, at the 
edge of such a deep gap, the gas density changes over a radial distance comparable to the disk scale height, while the 
width of the accretion band onto the circumplanetary disk is much narrower than the scale height (Fig. [3]). Therefore, 
we think that the gap formation would not affect the qualitative feature of the accretion flow. Another possible effect 
of gap formation is the disk size. As described in Ej4.ll several mechanisms are p ropos ed to explain the disk size. If 
the disk size is determined by the tidal effect as suggested by iMartin fc Lubowl (| 2 1 If ) . it should not be affected by 
gap formation, except when the Bondi radius is smaller than the Hill radius. On the other hand, if the disk size is 
determined by the radially outward velocity of the gas as described in §4.1[ lower density around the Hill sphere and 
thus stronger radial pressure gradient near the disk edge may enhance the outflow from the circumplanetary disk, and 
the disk size may become smaller. In any case, effects of gap formation should be examined in future works to check 
the validity of the results shown in this paper. 

5. SUMMARY 

In order to understand the structure of circumplanetary disks, we performed high-resolution hydrodynamic simula- 
tion and analyzed gas accretion flow in detail. We confirmed that gas accretion onto circumplanetary disks occurs in a 
manner that the gas is accreted from h igh altitude toward the disk surface do wnward with large incident angl e, which 
was suggested by previous studies (e.g.. lD'Angelo et al.] [2003; B ate et ai1l2003[ ). whereas lAvliffe fc Bate! (|2009bft showed 
that, in terms of mass flux, the downward accretion flow is not significant, which is inconsistent with our results. We 
found that the gas that has passed through the shock surface moves inward because its specific angular momentum is 
smaller than that of Keplerian rotation, whereas gas accretion through the midplane does not occur. While the net gas 
accretion across the Hill sphere is inward, outflow through near the two Lagrangian points is observed, although the 
regions around the two Lagrangian points have been, from the point of view of the potential energy, thought as a main 
accretion channel from protop lanetary disks. This outflow was not o bserved in previous hydrodynamic simulations 
for Jupiter-sized planets (e.g., iBate et all 120031 : lAvliffe fe~B atc 2009b). Outward radial velocity of the gas near the 
midplane significantly increases at a radial distance of about 0.2 times the Hill radius from the planet, and the gas 
can escape from the Hill sphere within a short period of time. 

We also obtained the distribution of mass and angular momentum of accreting gas onto the surface of circumplan- 
etary disks. We found that the accretion rates of mass and angular momentum can be well described by power-law 
functions. This distribution would be useful in the study of satellite formation, fo r ex ample, a radially one-dim ensional 
viscous-evolution model for circumplanetary disks, such as lWard fc Canupl (|2010[ ) and lMartin fc Lubowl (1201 1[). Recent 
development of viscou s modeling for MRI turbulence in protoplanetary and circumplanetary disks (jOkuzumi fc Hirosd 
l2011tlFuiii et al.ll2011[ ) would also contribute to construct more realistic models for the circumplanetary disks. However, 
in order to understand satellite formation processes, long term evolution of circumplanetary disks, which is determined 
by the evolution of accretion rate from protoplanetary disks to circumplanetary disks, is necessary and thus global 
evolution of protoplanctary disks with embedded giant planet s until the complete dissipation of protoplanet ary disks 
needs to be examined. Together with such global models (e.g.. lTanigawa fc Ikomall2~007t iWard fc Canuq|2010l ), we will 
be able to obtain long-term evolution of circumplanetary disks, which would provide better understanding of satellite 
formation processes. 

Our results demonstrate that gas accretion toward and within circumplanetary disks has a complicated vertical 
structure. The width of the accretion band toward a circumplanetary disk depends on the initial height of gas 
elements (Fig. [3]). Also, after accretion onto the disk, the gas just below the shock surface migrates inward, while the 
gas near the midplane moves radially outward (Fig. II T[) . Such vertical heterogeneity of the flow may have significant 
influence on the dynamical evolution of solid bodies in the circumplanetary disk, which we will examine in our future 
work. We also need to address the effect of viscosity and local calculation, which would affect the flow quantitatively. 
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Fig. 1. — Streamlines starting from four different heights (zo = 0.0,0.5,1.0,1.5), with xq = [±2,±3] and yo = ±L y /2. Interval of the 
starting points is 0.05 in the ^-direction. Green, blue, and red curves show streamlines in the horseshoe, accretion, and passing region, 
respectively. 
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Fig. 2. — Streamlines in the midplane of a circumplanctary disk. Left panel shows a streamline starting from R = 0.2 in the midplanc 
(£o,3/0,2o) = (0.2,0,0). Dashed lino shows contour line of <& = 0, which passes the two Lagrangian points Li (—1,0,0) and L2 (1,0,0). 
Colors show potential $; red regions are <E> > and blue regions are $ < 0. Right panel shows a closer view of three streamlines starting 
from R = 0.2 (same as the left panel), 0.1, and 0.05, respectively. 
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Fig. 3. — Classification of starting points of streamlines on the x—z plane at y = L v /2. Streamlines starting from the red region reach 
the boundary y = —L y /2 with x > (passing region), those in green region make U-turn and reach the boundary y = L y /2 with x < 
(horseshoe region), and those in blue region become trapped in f < 0.2 (accretion region). 
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Fig. 4. — Circumplanctary disk structure in the R-z plane with three different spatial scales. Log density is shown with colors, specific 
angular momentum is shown with contour lines, and gas velocity is expressed with arrows. All quantities are averaged in the azimuthal 
direction. 
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Fig. 5. — Mass flux on spheres with radii f = 1.0, 0.3, 0.1, and 0.03. Horizontal axis is azimuth angle and vertical axis is elevation 
angle on the sphere. Outward flux is defined to be positive. (<j>, 9) = (0, 0) and (180, 0) correspond to the sub-solar and anti-solar points, 
respectively. 
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Fig. 6. — Left panel shows azimuthally integrated mass flux at the radii f = 1.0 (solid red line), 0.3 (dashed green line), 0.1 (dot-dashed 
blue line), 0.03 (dotted purple line) as a function of elevation angle 9. Right panel shows azimuthally averaged radial velocity at the same 
radii. 
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Fig. 7. — Azimuthally averaged radial velocity in the midplanc normalized by local Keplerian velocity Vr, which corresponds to angle 
between the velocity vector of the flow and the circular orbit at the point. 
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Fig. 8. — Azimuthally averaged specific angular momentum at the midplane of the circumplanetary disk. Solid line shows specific angular 
momentum measured in the rotating frame of the Hill coordinate, dashed line is the case when it is measured in the inertial frame, and 
dotted line shows specific angular momentum of Keplerian rotation around the planet. 
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Fig. 9. — Azimuthally-averaged density at the midplane (left). Azimuthally averaged surface density (right), both as a function of the 
radial distance in the midplane. 
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Fig. 10. — Quantities along two streamlines starting from two different heights. Horizontal axis d is distance along each streamline from 
the starting point. Starting points for upper and lower panel are (2.37, L y /2, 0.0) and (2.37, Ly/2, 1.0), respectively. Thick solid (red) 
line shows Bernoulli integral, thick dashed (green) line shows kinetic energy, thick dotted (blue) line shows potential energy, and thick 
dot-dashed line (purple) shows logarithm of density. Thin dotted, dashed, dot-dashed lines show x, y, z, respectively. The streamline 
shown in the lower panel is accreted onto the circumplanctary disk, while that in the upper panel is not. Note that the streamline in the 
case of 2o = 1.0 is truncated at a point where the gas element starts rotating around the planet in the circumplanctary disk. 




Fig. 11. — Close up view of Fig. 1101 In addition to the quantities shown in Fig. 1101 azimuth angle <j> divided by it (thin dashed blue), 
specific angular momentum (thin dot-dot-dashed green), and logarithm of f (thin solid red) are also shown. 
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Fig. 13. — Azimuthally averaged mass fluxes through the surface z = z B (R) as a function of radius. Solid line shows / s and dashed line 
shows /Kep- Dotted line shows / s but excluding the region where / s < from the integration Eq. 1171 . f s is a flux directly observed at 
the surface, while /Kep is an effective flux assuming redistribution of radial location where specific angular momentum matches that of the 
local Kcplcrian rotation, (see Eqs. II17H and H21I 0. 
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Fig. 14. — Mass accretion rate onto the disk surface within a circle of radius. Solid line shows M B (Eq. H8\ ) and dashed line shows Mkcp 
(Eq. Dotted line shows M B but excluding the region where / s < from the integration of Eq. ]18[> . 
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Fig. 15. — Azimuthally-avcragcd specific angular momentum of the accretion flow onto circumplanctary disks denoted by j ZyB in the text 
(solid line). Dashed line is jKep = V ' wn i c h * s specific angular momentum for Kcplcrian rotation around the planet. 
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Fig. 16. — Mean specific angular momentum of the gas accreting onto the part of the circumplanctary disk within the radial distance R. 




Fig. 17. — Schematic picture of flow structure of circumplanctary disks. 



